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Abstract 

In this paper we present a new parallel algorithm for the solution of 

' ^ ' Hamilton-Jacobi-Bellman equations related to optimal control problems. 

'^S The main idea is to divide the domain of computation into subdomains 

following the dynamics of the control problem. This results in a rather 
T complex geometrical subdivision, but has the advantage that every sub- 

'"pi domain is invariant with respect to the optimal controlled vector field, so 

that we can compute the value function in each subdomain assigning the 
task to a processor and avoiding the classical transmission condition on 

' ' the boundaries of the subdomains. For this specific feature the subdo- 

mains are patches in the sense introduced by Ancona and Bressan in [l]. 
Several examples in dimension two and three illustrate the properties of 

Jr the new method. 
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^ 1 Introduction 

T— I 

^ The numerical solution of partial differential equations obtained by applying the 

Dynamic Programming Principle (DPP) to nonlinear optimal control problems 
is a challenging topic that can have a great impact in many areas, e.g. robotics. 
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aeronautics, electrical and aerospace engineering. Indeed, by means of the DPP 
one can characterize the value function of a fully-nonlinear control problem 
(including also state/control constraints) as the unique viscosity solution of a 
nonlinear Hamilton-Jacobi equation, and, even more important, from the so- 
lution of this equation one can derive the approximation of a feedback control 
for the system (see the next section for more details). This result is the main 
motivation of a PDE approach to control problems and represents the main 
advantage over other methods, such as those based on the Pontryagin minimum 
principle, that gives only an open-loop control (see e.g. [21] )■ It is worth to 
mention that the characterization via the Pontryagin principle gives only nec- 
essary conditions for the optimal trajectory (the state) and optimal open-loop 
control (the co-state). Although, from the numerical point of view, this sys- 
tem can be solved via shooting methods for the associated two point boundary 
value problem, in real applications the initial guess to start with is a particularly 
difficult choice for the co-state and often requires many efforts and tentatives. 
This is why it can be rather efficient to combine the dynamic programming and 
the Pontryagin approaches and use the approximate value function in order to 
compute a suitable initial guess for direct methods as proposed in [3T] . 

In this paper we mainly focus on the minimum time problem, which is asso- 
ciated to the following Hamilton- Jacobi-BcUman equation 

max{-/(x,a) • Vw(x) - 1} = 0, xeW'\no 

aeA (1) 

u{x) = , a; G ilo 

where d is the dimension of the state, A C K'" is a compact set defining the 
admissible controls, fig is the target set to be reached in minimal time and 
/ : R'' X yl — > M'' is the dynamics of the system. The value function m : M'' K at 
the point x is the minimal time to reach the target starting from x {u{x) — -f oo 
if the target is not reachable). For numerical purposes, the equation is solved in 
a bounded domain 17 D f^o- Boundary conditions on 9fl will be detailed later. 
Let us just mention that the same technique can be applied to more general 
control problems (see the last section for an example). 

The techniques used to obtain a numerical approximation of the solution of 
equation ([I]) have been mainly based on Finite Differences [HI [53] and Semi- 
Lagrangian schemes [221 [2S] • More recently, Finite Elements methods based on 
Discontinuous Galerkin approximations have been proposed, due to their ability 
to deal with non regular functions, which are the typical solutions appearing in 
the framework of viscosity solutions (see [T71I3S1 [H] ) . It is rather important to 
note that traditional approximation schemes presented for example in 22] and 
|18j are based on a fixed point iteration scheme, which computes the solution at 
each node of the grid at every iteration. Denoting by M the number of nodes in 
each dimension and considering that the number of iterations needed for conver- 
gence is of order 0(M), the total cost of this full-grid scheme is 0{J\1'^^^). We 
easily conclude that this algorithm is very expensive when the state dimension 
is d > 3, although it is rather efficient for low dimensional control problems as 
shown in [22] (see also the book |25j). 
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The "curse of dimensionality" issue has been attacked from several directions 
and new techniques have been proposed to accelerate convergence and/or to re- 
duce the memory allocation requirements. In [9j authors proposed an algorithm 
that allows to allocate only a small portion of the grid at every iteration. An- 
other proposal to reduce the computational effort is the so-called Fast Marching 
method [321 131] ■ While the full-size grid is always allocated, the computation is 
restricted to a small portion of the grid, thus saving CPU time. The cost of this 
method is of order 0{M'^ log M"^). In the original version, the Fast Marching 
method was derived for the Eikonal equation, corresponding (under a suitable 
change of variable) to equation ([T]) with f{x, a) = a and A = 5^(0, 1), the unit 
ball in E'' centred in (see [20j for details). Despite the fact that the Fast 
Marching method is very efficient, at present its application to more general 
equations of the form H{x,u{x),yu{x)) = is not an easy task and it is still 
under investigation (see [121 HSl [HI |35] ) . For these reasons other methods have 
been proposed exploiting the idea that by applying the same method in a finite 
number of pre-defined directions one can accelerate convergence. These "sweep- 
ing" methods have been shown to be efficient for the eikonal equation and, 
more recently, for rather general Hamiltonians [37j . 

A third possible strategy is based on the decomposition of the domain f2. 
The problem is actually solved in subdomains fi-', j = 1,...,R, whose size 
is chosen in order to reduce the number of grid nodes to a manageable size, 
see |26| . Therefore, rather than solving a huge problem in dimension d, one 
can solve R smaller subproblems working simultaneously on several processors. 
This produces a simple parallel algorithm. Depending on the choice of the 
subdomains il^ we can have some overlapping regions or a number of interfaces 
between the subdomains. The presence of interfaces and/or overlapping regions 
is a delicate point, since at each iteration of the algorithm it will be necessary 
to exchange information between processors to properly define the values on 
the interfaces. Without this communication load the result will not be correct. 
The interested reader can find in the book [32] a comprehensive introduction 
to domain decomposition techniques, whereas for an application to Hamilton- 
Jacobi equations we refer to [251 [12 • 

Finally, a decomposition of the domain based on the concept of "patchy 
feedbacks" has been proposed. The name "patchy" has been introduced in a 
work by Ancona and Bressan [1], where the authors studied the problem of the 



asymptotic stabilization of a control system. Their main result (see Theorem 2.1 



in Section 2.3 1 states that, under suitable assumptions on the control system, 
stabilization can be obtained by means of a special feedback control, the patchy 
feedback, which is piecewise constant on a particular partition of the domain. 
Such a partition has the fundamental property that each part, or "patch", is 
invariant with respect to the optimal dynamics driving the system. This re- 
sult can be exploited from a numerical point of view. Indeed, once a patchy 
decomposition has been obtained, the computation can take advantage of the 
invariance of the subdomains assigning each patch to a processor, yielding a very 
efficient parallel algorithm that drops the need of communication between the 
processors. This is the spirit of what we call a "patchy method" . Unfortunately, 
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the result of Ancona and Bressan is purely theoretical and their patchy decom- 
position turns out to be not constructive. Then, one has to face the problem of 
a numerical approximation of a dynamics invariant domain decomposition. 

A first example of discrete patchy method has been proposed by Navasca and 
Krener in [301 |3T]. The authors adopt a formal method developed by Al'brekht 
[3] that essentially translates the Hamilton- Jacobi-Bellman equation associated 
to a control problem into a system of algebraic equations, whose unknowns are 
the coefficients of the developments in power series of the cost and optimal 
feedback of the control problem. This gives an approximate solution in a small 
neighbourhood of the origin, which is the first patch of their domain decom- 
position. The solution is then extended in new patches around the first, by 
picking some boundary points from where optimal trajectories emanate (com- 
puted numerically backward in time). Those points define the centers of new 
neighbourhoods that can be used to restart the method. The solution is then 
obtained iteratively by fitting together the approximations in all the patches. 
More recently, it has been shown that this technique can be extended to obtain 
high-order accuracy in the regions where the value function is smooth (see [2 8) 
for more details). 

Despite the high speed of the method, that actually does not use any grid, 
there are many open questions on its application. The first limitation is the 
strong regularity assumptions on the solutions necessary to set the problem in 
these terms. Indeed, it is well known that the most simple control problems may 
have optimal feedbacks which are not even continuous (see [27] and [S]). The 
second crucial point is the construction of the patchy decomposition that, in the 
examples contained in [30|, ISlj , does not result in the invariance with respect 
to the optimal dynamics. This makes the solution rather inaccurate, especially 
near the boundaries of the patches. 

The goal of this paper is to present and test a new patchy technique based 
on the coupling between a semi-Lagrangian scheme and a domain decomposi- 
tion leading to a dynamic partition of the domain into subdomains which 
are invariant with respect to the optimal dynamics and to the optimal trajecto- 
ries to the target. To this end we will use the feedback controls computed by 
means of the semi-Lagrangian scheme on a rather coarse grid, then this informa- 
tion will be plugged into a dynamic algorithm which allows to obtain a domain 
decomposition where the domain boundaries correspond to approximate opti- 
mal trajectories. Finally, the computation of the value function on a fine grid is 
obtained via a parallel algorithm which assigns every sub-domain to a processor. 

The paper is organized as follows. Section 2 is devoted to the general presen- 
tation of our problem, of the semi-Lagrangian scheme and of classical domain 
decomposition techniques for Hamilton- Jacobi equations. Moreover, we briefly 
describe the results on patchy methods which have been proved by Ancona 
and Bressan [T]. In Section 3 we present the patchy domain decomposition 
method and our algorithm to split the domain into invariant subdomains. We 
discuss there several issues related to the implementation of the method and its 
parallelization. Finally, Section 4 is devoted to the numerical tests on control 
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problems in dimension 2 and 3. 



2 Background 

In this section we briefly introduce the numerical scheme we adopt to discretize 
equation ([l]) and the classical domain decomposition technique for Hamilton- 
Jacobi equations, including an algorithm that will be used in the following. 
Next, we recall the notion of patchy decomposition and the result by Ancona 
and Bressan concerning the asymptotic stabilization of control systems by means 
of patchy feedbacks, that inspired our patchy numerical method. We refer the 
interested reader to [H [21 |3j for details. 

2.1 The semi-Lagrangian scheme 

Let us denote by > a fictitious time step (see the book [5S] for details) and 
by fc > the space step. We introduce a structured grid G on il with nodes 
Xi, i = 1, . . . , N. We also denote by G the internal nodes of G and by dG its 
boundary, whose nodes act as ghost nodes. We map all the values at the nodes 
onto a A'^-dimensional vector U = {Ui, . . . , Un)- By a standard semi-Lagrangian 
discretization [6l [Tj [22] of ([T]), it is possible to obtain the following scheme in 
fixed point form 

U = F{U) , (2) 
where F : R^ is defined componentwise by 

f min {I [U] {x, + hf{x„ a)) + h} x,eG\no, 

[FiU)]^=< x.eQo, 
[ +00 Xi € dG . 

The discrete value function U is extended on the whole space f2 by a linear 
d-dimensional interpolation, represented by the operator I as described in |14j . 
We choose a variable time step h = hi^a such that \hi^af{xi,a)\ ~ k for every 
i = 1, ... ,7V and a € A, so that the point Xi + hi^af{xi,a) falls in one of the 
first neighboring cells. The minimum over A is evaluated by direct comparison, 
discretizing the set A with Nc points. Note that the definition F{U) — +oo 
on dG corresponds to impose state constraint boundary conditions. The final 
iterative scheme reads 

We refer to [551 US] for details and convergence results. It is important to note 
that the scheme (|2| produces an approximate feedback at every node of the 
grid. Moreover, the knowledge of the approximate value function at the nodes 
of the grid allows to extend it everywhere in f2 via interpolation. Then we can 
always obtain a feedback map '■ ^ ^ A just defining 

^h{x) :=argmin{I[C/](a;, + /i/(xi,a)) + /i} (4) 
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Under rather general assumptions (see [53]), it can be shown that this is an 
approximation of the feedback map constructed for the continuous problem as 

$(a;) :— argmax{— /(a;, a) • \7u{x) — 1} (5) 

A detailed discussion on the construction of feedback maps via the value function 
is contained in jS] p. 140-143. It is important to note that weak convergence 
results apply also for Lipschitz continuous value functions. 



2.2 Domain Decomposition method 

The domain decomposition method allows to split the problem in fl into R 
subproblems in subsets $7-', j = 1, . . . ,R such that fl — il^ U . . . U il^. For 
every pair j ^ I oi indexes corresponding to adjacent subdomains U,^ and fi^, 
let us denote by fi-'^ the non-empty overlapping zone il^ H which is assumed 
to contain at least one grid cell. This decomposition of the domain induces 
a decomposition of the N degrees of freedom, so that to each subdomain Q,^ 
it is associated a number of nodes This allows to consider the restriction 
of U to U,^ that we denote by with j — 1, . . . ,R. Similarly, we can define 
in a natural way the restriction to the subdomain W of the operator in ([2]), 
denoted by : R^' -> M^' . Then, we can define globally in the following 
splitting operator Fsplit : — > K^, given componentwise by 

r [F^iW)]i if 3j such that x, e \ U^^j , 

[i^sPLiT(?/)]i = < {[F^W)],} otherwise. 

(6) 

Following dS] it is easy to prove that fixed point iterations for F and -Fsplit 
have the same solution. 

We now describe a simple algorithm to compute the fixed point of Fsplit- 

Domain Decomposition Algorithm: 

Step 1. (Initialization) For n — the initial guess ;7(°) e R^ is fixed to on the 
nodes corresponding to the target fig and +oo elsewhere. 

Step 2. (Computation) is computed separately and in parallel in every 

subdomain D,^ by 

f^j,(«+l/2) ^ pj^ljj.in)^ j = 1, . . . 



Step 3. (Coupling and synchronization among processors) 



C/,f'^"+^/^' if 3j such that X, e \ Uf^j , 

min I 1 otherwise . 

(7) 
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step 4. (Stopping criterion) If |lC/("+^)-t/(")|loo >toll go to Step 2 with n ^ n+1, 
otherwise stop. 

In order to speed up the convergence of the algorithm above, we use iterations 
of Gauss-Seidel type, meaning that we employ the updated values of the nodes 
as soon as they are available. From now on the final solution computed by the 
Domain Decomposition algorithm will be denoted by fee. 

2.3 Patchy feedbacks and stabilization for controlled sys- 
tems 

The concept at the basis of the numerical method we present in the next section 
is the notion of patch. It has been introduced for the first time by F. Ancona 
and A. Bressan [I], [2] in the context of the stabilization of controlled systems. 
Let us recall here for completeness their main definitions and results. 

We consider the control system 

m = fiyit),ait)) ait)eA, (8) 

assuming that the control set A C is compact and the dynamics f •.'R'^xA ~> 
M.'^ is sufficiently smooth. Moreover we choose as admissible controls all the 
functions a(-) belonging to 

A := {a : (0, oo) — > A : a(-) is measurable} . 

For every initial point x € M'^ and admissible control ao E A, we denote by 
2/(- ; X, ao) the corresponding trajectory, which is an absolutely continuous func- 
tion defined on some maximal interval [0; t"'^'^^ (y)) , satisfying the system ([s]) for 
a.e. t > with initial condition y{0) ~ x and control cq. 

The following definition extends to control systems the classical notion of sta- 
bility. 

Definition 2.1 The system ^ is said to he globally asymptotically controllable 
(to the origin) if the following holds: 

1. for each x G M'' there exists some admissible control Oq such that the 
trajectory t y{t) — y{t;x,ao) is defined for all t > and y{t) — >■ as 
t — > oo, 

2. for each £ > there exists ^ > such that for each x € M'^ with \x\ < 5 
there is an admissible control Oq as in 1. such that \y{t)\ < e for allt>0. 

Given an asymptotic controllable system, a classical problem is to find a feed- 
back control a — k{y) : — > A such that all the trajectories of the correspond- 
ing closed loop system 

v^f{y,k{y)) (9) 
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tend asymptotically to the origin. 

Since this problem may not admit any solution in the class of continuous feed- 
backs, Ancona and Bressan introduce and investigate the properties of a partic- 
ular class of discontinuous feedbacks, the so-called patchy feedbacks. 

The following definition gives the fundamental concept of patch. 

Definition 2.2 Let be fl C M.'^ an open domain with smooth boundary dQ and 
g a smooth vector field defined on a neighborhood of We say that the pair 
{^l,g) is a patch if Q, is a positive-invariant region for g, i.e. at every boundary 
point y £ dfl the inner product of g with the outer normal n satisfies 

{9{y).n{y)) < 0. 

Then, by means of a superposition of patches, we get the notion of patchy vector 
field on a domain il C M''. 

Definition 2.1 We say that g : Q, ^ is a patchy vector field if there exists 
a family of patches {(ri",^") : a G X} such that 

- X is a totally ordered index set, 

- the open sets form a locally finite covering of fi, 

- the vector field g can be written in the form 

5(y) = 5"(y) y£n^\[jnP. 

We use (fl,g, , g°')aei) to indicate the patchy vector field and the family of 
patches. 

By applying the previous definitions to the closed loop system ^ we define a 
patchy feedback control as a piecewise constant map k : M.'^ — > A such that the 
vector field g{y) := f{y, k{y)) is a patchy vector field. More precisely: 

Definition 2.2 Let (fl, g, (Q" , g°')a^x) be a patchy vector field. Assume that 
there exist control values k" £ A such that, for each a £ X 

/3>Q 

Then the piecewise constant map 

k{y) = A" if y e r!" \ y f]'' 

I3>a 

is called a patchy feedback control on VL. 
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Definition 2.3 A patchy feedback control k : M'^\{0} A is said to asymptoti- 
cally stabilize the closed loop system ^ with respect to the origin if the following 
holds: 

1. for each x G M"^ \ {0} and for every trajectory y{-) of ^ starting from x 
one has y{t) as t ^ T^^^{y), 

2. for each e > there exists (5 > such that, for each x G \ {0} 
with \x\ < S and for every trajectory y{-) of ^ starting from x one has 
\y{t)\ < e, for allO<t< r'^'^'^iy). 

Finally, the main result of Ancona and Bressan can be summarized as follows. 

Theorem 2.1 // the system ^ is asymptotically controllable, then it admits 
an asymptotically stabilizing patchy feedback control. 



3 The patchy domain decomposition 

In this section we introduce our new numerical method for solving equations 
of Hamilton-Jacobi-Bellman type. In particular we focus on the minimum time 
problem ([T]). The main feature of the new method is the technique we use to 
construct the subdomains of the decomposition, which are "patches" in a sense 
inspired by the definitions of the previous section. Indeed, we will see that these 
patches turn out to be quite invariant with respect to the optimal dynamics 
driving the system. Moreover, their boundaries can be rather complicated, but 
this is not a major difficult for the technique, since we do not need to apply any 
transmission condition between them. 

Let us introduce two rectangular (structured) grids. The first grid should be 
rather coarse because it is used for preliminary (and fast) computations only. It 
will be denoted by G and its nodes by xi, . . . , Jj^i where N is the total number 
of nodes. We will denote the space step for this grid by k and the approximate 
solution of the equation ([l]) on this grid by Up. 

The second grid is instead fine, being the grid where we actually want to 
compute the numerical solution of the equation. It will be denoted by G and 
its nodes by Xi,. . . ,xn, where N is the total number of nodes {N»N). We 
will denote the space step for this grid by k and the solution of the equation 
([!]) on this grid by Up. We also choose the number R of subdomains (patches) 
to be used in the patchy decomposition and we divide the target Oq in R parts 
denoted by ^q, with j = 1, . . . , R. 
The patchy method can be described as follows. 

Patchy Algorithm: 
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Step 1. (Computation on G). We solve the equation on G by means of the classical 

For coherence 



domain decomposition algorithm described in Section 2.1 



we choose the (static) decomposition made by R subdomains (as the num- 
ber of patches). This leads to the function Up. 



Step 2. (Interpolation on G). We define the function Up on the fine grid G 
by interpolation of the values Up. Then, we compute the approximate 
optimal control 

Z*{xi) = a.Tgm.ln{l[U'^°\x, + h,^af{x^,a)) + h,^a}, x,&G. (10) 

Even if a* is defined on G, we still use the symbol "tilde" to stress that 
optimal controls are computed using only coarse information. We delete G. 



Step 3. (Main cycle) For every j = 1, . . . , R, 



Step 3.1. (Creation of j-th patch). Using the (coarse) optimal control a*, we 
find the nodes of the grid G that have the part JIq of the target in 
their numerical domain of dependence. This procedure defines the 
j-th patch, naturally following the (approximate) optimal dynamics. 
This step will be detailed later in this section. 

Step 3.2. (Computation in j-th patch). We initialize the j-th solution equal 
to +00 on the j-th patch and equal to on the part ftf^ of the tar- 
get (initial guess). Then, we apply iteratively the scheme ^ in the 
j-th patch until convergence is reached. Finally, the j-th solution is 
copied in the matrix that will contain the global solution Up. 



Details on Step 3.1. Following TTl, the basic idea we adopt here is to divide 
the whole domain starting from a partition of the target only, letting the dy- 
namics compute the partition in the rest of the domain. To this end we use 
the approximation of the optimal control given by a* and then we obtain a 
domain decomposition fully compliant to the dynamics. More precisely, we di- 
vide the target Q.q in R parts, each associated to a colour indexed by a number 
j = 1, . . . , i?. Assume for instance that Oq is a ball at the centre of the domain 
and focus on the subset of the target with a generic colour j, denoted by I^q, 
see Fig. ^a). The goal is to find the subset of the domain which has JIq as 
numerical domain of dependence. To do that, we initialize the grid nodes with 
the values (pi as follows 
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Then we solve the following ad hoc scheme, similar to the fixed-point algorithm 
([2| for the main equation 

(^i = l[(l)]{xi + h,f{xi,a*{xi))) , x,&G. (11) 

Here hi is chosen in such a way that \hif{xi,'d* {xi))\ — k. Once the computation 
is completed, the whole domain will be divided in three zones: 

A{ = {xi : = 1} , = {x, : = 0} , A| = {x, : (t>i £ (0, 1)} , 

see Fig. [ijb). This is due to the fact that the interpolation operator I in the 



semi-Lagrangian scheme ( 11 ) mixes the values (pi through a convex combination, 
thus producing values in [0, 1] even if the initial datum is in {0, 1}. Since we 
need a sharp (non-overlapping) division of the domain, we "project" the colour 
j into a binary value 



1 , > i 

, (|)^<l 



X, e G. (12) 



and then we define the subdomain = {xi £ G\Qi : 4>\ = 1} as the j-th 
patch, see Fig. [ijc). Once all the patches j — l,...,i? are computed, they 
are assembled together on the grid G. Thus the grid results to be divided in R 
patches, each associated to a different colour, as shown in Fig. [ijd). 

The main point here is that patches 's are constructed to be invariant with 
respect to the optimal dynamics, meaning that the solution of the equation in 
each patch will not depend on the solution in other patches. This is equivalent 
to state that there is no crossing information through the boundaries of the 
patches. 

We stress that Step 3.1 of the algorithm is not expensive, even if it is per- 
formed on the fine grid G. The reason for that is the employment of the pre- 



computed optimal control a* in the equation (11), which avoids the evaluation 



of the minimum (see the scheme ([2j)). Moreover, the stop criterion for the fixed 



point iterations (11) can be very rough, since we project the colors at the end 



and then we do not need precise values. 

Remark 3.1 (boundary conditions). To make it evident that each patch 
is independent of the others, we impose state constraint boundary conditions 
on the boundary of the patches. This kind of boundary conditions force the 
optimal direction /(x, a*) to point inside the patch. If patches are not invariant 
with respect to the optimal dynamics, this boundary condition leads to a huge 
error, which propagates all over the domain. It is important to note that, if the 
optimal trajectory runs along the boundary of a patch and the boundary of the 
patch is not aligned to the grid, the semi-Lagrangian scheme does not allow to 
select an optimal control which drives the dynamics exactly along the boundary, 
even if the set of admissible directions /(x. A) include it. This is caused by the 
interpolation operator which necessarily makes use of some nodes outside the 
patch, where state constraints are found. The result is selecting an optimal 
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(a) 





(d) 

Figure 1: Creation of patches for a test dynamics, R—A, rio^small ball in the 
centre: (a) Select a subdomain fl^ of the target ^Iq. (b) Find the nodes which 
depend, at least partially, on fig. (c) Define fl^ projecting the color in a binary 
value, (d) Assemble all patches. 



direction which points toward a cell fully inside the patch (and not across two 
patches), see Fig. [3]-(d) and its caption. 

Remark 3.2 (patches as a partition ofG). We have no guarantee that patches 
ri-^'s do not overlap nor that they cover the whole domain. Over the overlapping 
zones we can simply choose a colour at random. Instead, if they do not cover all 
the domain we can repeat the computation in the not-coloured nodes relaxing 



the condition in ( 12 ), i.e. choosing a different value for 1/2. Alternatively, in the 
case of isolated not-coloured points, we can assign to them the colour of their 
neighbours. 

Remark 3.3 (equivalence of patches) . Patches are not meant to be equivalent 
in terms of number of nodes, nor in term of CPU time needed to compute the 
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solution in them. As we will see in the following, the difference in CPU time 
between the patches will be the key property that lets the patchy algorithm 
overcome the classical domain decomposition method. On the other hand, too 
large differences will lead to a deterioration of the performance. 

Strategies for parallelization 

The patchy algorithm can be parallelized in two ways. 

Method 1 : Patches are processed one after the other and only the compu- 
tation in each patch is performed in parallel, assigning a batch of nodes among 
processors. This strategy gives priority to saving CPU time. Indeed, processors 
are active all the time. On the other hand, the full grid G must be allocated in 
RAM and be visible by all processors all the time. 

Method 2: Patches are distributed among processors and computation of 
each patch is serial. This strategy gives priority to memory allocation issues. 
As patches are invariant with respect to the optimal dynamics, processors do 
not need to communicate and they can have a non-shared memory. Each pro- 
cessor works on its patch and once the computation is done it returns the result 
to the master thread. Note that if the number of processors equals the number 
of patches this procedure is not efficient in term of CPU time, because once a 
processor finished its job, it can not be assigned to another one and remains 
idle. Instead, if the number of processors is less than the number of patches, 
processors are better used because once one of them has finished its job, it can 
take another job until all the patches are covered. 

Remark. All tests presented in this paper are performed implementing 
Method 1 and in the following we will always refer to this choice. 

4 Numerical investigation in dimension two 

In this section we first list the dynamics considered for the numerical tests. 
Then, we investigate the optimality of the patchy decomposition and the per- 
formance of the algorithm with respect to the classical domain decomposition. 

Numerical tests were performed on a server Supermicro 8045C-3RB using 1 
CPU Intel Xeon Quad-Core E7330 2.4 Ghz and 8 x 4 GB RAM running under 
Linux Gentoo operative system. 

4.1 Choice of benchmarks 

We will test the method described above against three minimum time problems 
of the form ([l]). The numerical domain is always f2 — [—2, 2]^. 



13 



Test 1 (Eikonal) 



d=2, f{xi,X2,a)^a, ^ = ^2(0,1), ^io = ^2(0, 0.5). 
Test 2 (Fan) 

d = 2, /(a;i,X2,a) = |xi+a:2+0.1|a, A = 32(0,1), f^o = {a^i = 0}. 
Test 3 (Zermelo) 

d=2, f{xi,X2,a) = 2.1a + {2,0), ^ = ^2(0,1), r!o - -82(0, 0.5). 

In Fig. [2] we show the patchy decomposition for the three dynamics described 
above in the case R = 8, Nc = 32, N = 50 and N = 100. We also superim- 
pose the optimal vector field f{x,a*) to show that patches are indeed (almost) 
invariant with respect to the optimal dynamics. Only a few arrows cross from 
a patch to another. We note that patches cover the whole domain but they are 
not equivalent in terms of area, even if the target fio was divided in i? = 8 equal 
parts to generate the decomposition. 

4.2 Optimality of the patchy decomposition 

Beside the discretization error due to the numerical algorithm, the patchy algo- 
rithm introduces another error, which will be denoted by Ep. Error Ep can be 
found comparing the solution of the patchy algorithm with that of the classical 
domain decomposition method based on the same numerical scheme, 

Ep := Up — Udd- 

This additional error is exclusively due to the fact that patches are not com- 
pletely dynamics- invariant, because they are found using information computed 
on the coarse grid G. It is plain that we could in principle compute an optimal 
patchy decomposition for the grid G working directly on the fine grid G, but 
this would cost as much as computing the solution U on G. 

In the following we study the norms ||i?p||i and ||i?p||oo as a function of the 
space step size k and k. We report the results for R = 16, which is the largest 
number of patches we tested. Note that error Ep necessarily increases as R 
increases because the number of boundaries increases. We recall that we impose 
state constraints boundary conditions at the boundary of each patch, fixing a 
large value of the solution outside the patch. In this way each computation is 
completely independent of the others and the quality of the final global solution 
Up at the boundaries of the patches only depends on the degree of invariance 
of the patches with respect to the optimal dynamics. 

Results for the three dynamics are shown in Tables [T] [2] and [Sj 
We see that the first line of each table reports in many cases unsatisfactory 
results, caused by the excessive roughness of the grid G (see the case fc=0.08, 
corresponding to A^=50). Even the case k = k (i.e. the grid is not refined at all) 
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Figure 2: Patchy decompositions with i? = 8, TV^ = 32, = 50 and N = 100. 
Only few arrows are shown for clarity of visualization, (a) Eikonal, (b) Fan, (c) 
Zermelo, (d) a detail of (b). 



is not satisfactory. This can be explained by recalling that the computations 
on the two grids are not equivalent because the second one employs the state 
constraint boundary conditions at the boundaries of each patch. If the grid is 
not fine enough, the error due to the boundary conditions is large, and tends 
to propagate inside each patch. Conversely, if G has at least 100 nodes per 
dimension {k >0.04), the behaviour of the error is surprisingly good because it 
decreases as k decreases (for any fixed k) and ||-Ep||i is of the same order of k 
itself. Two comments about this behaviour are in order. First, the error caused 
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Table 1: Patchy error ||£;p||i (||i;p||co)- Dynamics Eikonal, = 32, R= 16 





k = 0.08 


k = 0.04 


k = 0.02 


k = 0.01 


k = 0.005 


K=0.08 


0.436 (0.960) 


0.275 (1.856) 


0.102 (0.048) 


0.065 (0.034) 


0.048 (0.026) 


^=0.04 




0.088 (0.046) 


0.029 (0.023) 


0.014 (0.042) 


0.005 (0.008) 


K=0.02 






0.038 (0.029) 


0.012 (0.013) 


0.004 (0.008) 


^=0.01 








0.011 (0.016) 


0.006 (0.010) 


^=0.005 










0.004 (0.008) 



Table 2: Patchy error (||£;p||oo). Dynamics: Fan, Nc = 32, R = 16 





k = 0.08 


k = 0.04 


k = 0.02 


k = 0.01 


k = 0.005 


K=0.08 


1.393 (3.023) 


0.123 (1.507) 


0.037 (0.315) 


0.017 (0.263) 


0.011 (0.263) 


K=0.04 




0.114 (1.502) 


0.032 (0.149) 


0.011 (0.095) 


0.006 (0.095) 


K=0.02 






0.032 (0.111) 


0.011 (0.061) 


0.004 (0.037) 


n=o.oi 








0.011 (0.079) 


0.004 (0.037) 


^=0.005 










0.004 (0.037) 



Table 3: Patchy error ||i^p||i (HE'pIIcx))- Dynamics: Zermelo, Nc = 32, R = 16 





k = 0.08 


k = 0.04 


k = 0.02 


k = 0.01 


k = 0.005 


^=0.08 


0.171 (0.293) 


0.159 (0.059) 


0.097 (0.057) 


0.026 (0.027) 


0.006 (0.016) 


^=0.04 




0.101 (0.063) 


0.033 (0.041) 


0.011 (0.023) 


0.004 (0.016) 


K=0.02 






0.039 (0.039) 


0.012 (0.023) 


0.004 (0.016) 


K=0.01 








0.011 (0.020) 


0.005 (0.015) 


^=0.005 










0.004 (0.016) 



by the state constraints boundary conditions stays close to the boundaries and 
does not propagate in the interior. Second, the error in locahzing the patches 
is much less effective than the error caused by the numerical scheme. Let us 
investigate the latter point in more detail: If patches are not at all invariant with 
respect to the dynamics, as in the classical domain decomposition algorithm, 
we do not expect the error decreases as k decreases, because in this case the 
error is mainly due to the missing information from the boundaries rather than 
to the numerical scheme. If instead patches are perfectly independent of each 
other, the error decreases as k decreases, because the error is now only due to 
the numerical scheme. Our numerical tests show a behaviour much more similar 
to the second case, in which the error of the numerical scheme is leading with 
respect to that of the decomposition. 

We also note that the L°° error is always larger than the error. Quite 
often we find a very small mimbcr of nodes with a large error near the boundaries 
of the patches, especially at those nodes where two patches and the target meet. 
This mainly affect the L°° error but not the error. Finally we note that the 
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results are similar for the three dynamics, showing a good robustness even for 
highly rotating vector fields like that of Fan dynamics. 

Fig. [3]-(a,b,c) shows the functions Ep for the three tests. These results 
confirm that the error is concentrated along the boundaries of the patches. In 
the Eikonal and Zermelo case the error starts propagating from the target and 
increases as long as characteristics go away. In the Fan case, instead, the largest 
error is found where patches and target meet. Note that in the Eikonal case (a) 
no error is found where patches boundaries are aligned to the grid, since the 
interpolation makes no use of nodes belonging to different patches. Fig. [3]-(d) 
shows a detail of the Fan decomposition along with the approximate optimal 
vector field computed by means of the final solution Up. The effect of the state 
constraints is perfectly visible (arrows point unnaturally inward) confirming the 
fact that each patch is computed independently (see Remark 3.1). 

Fig.|4]shows the value function U and its level sets for the Eikonal test. Here 
we see small perturbations where patches meet. It is interesting to note that we 
they not meet forming a discontinuity, but rather a hollow. 

4.3 Comparison of CPU times 

In this section we compare the patchy algorithm with the domain decomposition 
algorithm in terms of CPU time. Before that, we report in Table [3] the times 
(in seconds) for the single steps of the patchy algorithm. Times for the Step 3.1 
are the most interesting ones because this step is expected to be the slowest one 
after Step 3.2 (main computation on the fine grid). Thus, time spent in Step 
3.1 can completely neutralize the advantage we hope to get in the subsequent 
main computation. As we can see, Step 3.1 is much more costly than Steps 1 
and 2, but not so much compared with the main computation. 

Tables [Sj |6] [t] and [s] report the total CPU time (in seconds) for the three 



Table 4: Processors: 4, 7Vc=32, Grid: N = 100^ N = 800^, i?=16 





Step 1 


Step 2 


Step 3.1 (all j's) 


Step 3.2 (all fs) 


Eikonal 


2 


1 


23 


409 


Fan 


2 


2 


52 


796 


Zermelo 


2 


1 


30 


512 



dynamics of Sect. 4.1 as a function of the number of cores (1-4) and the number 
of patches (i?=2,4,8,16). For the Eikonal test we also vary the number of discrete 
controls {Nc=l& and 32). These results are compared with the best outcome 
of the domain decomposition method obtained varying the number of domains 
(again 2,4,8,160 

^The CPU time of the domain decomposition method does not vary a lot varying the 
number of domains, but small differences are present. They are due to the different order in 
which nodes are visited and synchronization overhead at the end of each iteration. 



17 



(a) 



(c) 



II 



(b) 



■ 
















































































































































































































































































































































/ 








































V 














/ 
















































/ 




































V 


























































/ 












































































































































































































/ 












r 














































/ 






































































































/- 










/ 










































V 








































































'A 
























































'- 






















'- 
















































































































A, 
















































































































t 






























































































i 






t- 


















t- 


r- 


! 










4- 













(d) 



Figure 3: Patchy error Ep, N = hi] ^ N = 100 for (a) Eikonal, (b) Fan, (c) 
Zermelo. In (d) it is shown a detail of the patchy decomposition for the Fan 
dynamics, together with the optimal vector field f{x,a*) computed by means 
of the final patchy solution. 



Table 5: Dynamics: Eikonal. Nc=l6. Grid: N = 100^ N = 800^ 





2 domains 


4 domains 


8 domains 


16 domains 


Best DD 


1 core 


1547 


1076 


1058 


933 


1571 


2 cores 


845 


595 


574 


504 


820 


4 cores 


459 


325 


317 


271 


415 
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(a,j (b) 

Figure 4: Patchy solution for the dynamics Eikonal. (a) Value function and (b) 
a detail of its level sets. Small hollows are visible in correspondence of the lines 
{x — y} and {x = —y} as in Fig. |3[(a). 



Table 6: Dynamics: Eikonal. Nc^32. Grid: N = 100^ N ^ 800^ 





2 domains 


4 domains 


8 domains 


16 domains 


Best DD 


1 core 


2702 


1897 


1843 


1623 


2785 


2 cores 


1462 


998 


968 


872 


1430 


4 cores 


771 


532 


514 


435 


716 



First, we see that the speed-up with respect to the number of cores is very 
satisfactory and proves that the parallelization Method 1 (see Sect. [3]) we im- 
plement here is sound. Second, we see that the CPU time decreases remarkably 
as the number of patches R increases. For i?=16 the CPU time is largely less 
than that of the best domain decomposition method. This is one of the main 
results of the paper. 

Differences among iVc=16 and iVc=32 are instead less clear, although the 
patchy algorithm should have an advantage for large Nc because of the smaller 
ratio between CPU time for Step 3.1 (one discrete control) and Step 3.2 {Nc 
discrete controls). 

In order to understand why patchy algorithm is faster than domain decom- 
position method, let us consider again the Eikonal case with R—8, see Fig.[2]-(a). 
If we visit the nodes in a single predefined order (i.e. we do not implement the 
Fast Sweeping technique [391 or similar ones), the eight subdomains need a dif- 
ferent number of iterations to reach convergence. This is due to the fact that 
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Table 7: Dynamics: Fan. iVc=32. Grid: N = 100^ -> iV = 800^ 





2 domains 


4 domains 


8 domains 


16 domains 


Best DD 


1 core 


3712 


3322 


3049 


3172 


4163 


2 cores 


2020 


1746 


1596 


1559 


2124 


4 cores 


1032 


900 


841 


852 


1069 



Table 8: Dynamics: Zermelo. iVc=32. Grid: N ^ 100^ N ^ 800^ 





2 domains 


4 domains 


8 domains 


16 domains 


Best DD 


1 core 


3113 


2675 


2126 


2018 


3209 


2 cores 


1651 


1404 


1111 


1054 


1640 


4 cores 


871 


721 


584 


545 


825 



for some of them the visiting order corresponds to the upwind direction, while 
for the other domains the visiting order corresponds to the downwind direction. 
If we do not know a priori that the eight domains are invariant with respect to 
the optimal dynamics, we can not label as "done" the computation in a domain 
before computations in all domains are fully completed, because in any moment 
a new information can enter the domain, making necessary new computations. 
On the contrary, if we know a priori that domains do not depend on each other, 
we can label as "done" the computation in a domain as soon as the solution 
reached convergence, and use computational resources for other domains. Note 
that the more the number of domains the more efficient is the computation. 

Finally note that usage of the Fast Sweeping technique can mitigate the per- 
formances of patchy method, but not neutralize them completely. Moreover, the 
patchy algorithm has the clear advantage that no synchronization nor crossing 
information among processors is needed. This is a great advantage when using 
distributed memory parallel computers, where communications are performed 
via cables connecting cluster nodes. This advantage is not really included in 
our experiments because our cores share a common RAM. 

4.4 Patchy method with obstacles 

We have also tried to run the patchy algorithm in a minimum time problem 
(Eikonal dynamics) with obstacles. In Fig. [5] we show the obstacles (the black 
circle and the rectangle) , the level sets of the solution, the patchy decomposition 
and the patchy error Ep. The behaviour of the patchy decomposition is good 
because the dynamics drives the patches around the obstacles. If not influenced 
by the obstacles, the error is concentrated around the boundaries of the patches 
as expected. Instead, when a boundary meets an obstacles, the error can either 
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stop propagating (see the circle) or spread out (see the rectangle, right side). 




(a) (b) 

Figure 5: Patchy domains bypass the obstacles driven by the dynamics, (a) 
Level sets of the value function and patchy decomposition, (b) Error Ep. 



4.5 Limitations of patchy method 

The overall efficiency of the patchy method depends on the dynamics and the 
shape of the target ilo. Unfortunately it is not always possible to reach a 
suitable patchy decomposition which allows to run the algorithm. This can 
happen for example if the target is very small and then it can not be divided in 
R subdomains. 

Another issue, much more difficult to fix, comes out whenever there is a large 
difference between the sizes of the patches and possibly some of them degenerate 
in a subset of a few grid nodes. This is the case of the classical "Lunar Landing" 
problem 

d = 2, f{xi,X2,a) = {x2,a), ^ = {-1,1}, no = B2(0,e). 

In this case the patchy decomposition consists of 2 large domains and R — 2 
smaller domains, see Fig. [6]-(a). The small domains degenerate to empty sets 
when e tends to zero, because all the optimal trajectories tend to meet in only 
two switching lines. 

A third dangerous case arises when some regions in fig are not reachable. 
If the dynamics make it impossible to reach the target from some point x € fi, 
the value function is set to u{x) — +oo. From the numerical point of view, 
the solution stays frozen at the value given as initial guess. At these points the 



optimal control ( 10 ) is not uniquely defined, and then the patchy decomposition 
can not be build. On the other hand, in the not-reachable regions the solution 
is in some sense already computed, then the issue can be easily fixed by a slight 
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modification of the algorithm. After Step 1 we locate the regions where Up is 
very large and then we do not consider those regions in the rest of computations. 



4.6 Patchy decomposition for non-target problems 

In the case of non-target problems we can not in principle build the patchy 
decomposition. Indeed, we recall that our patchy decomposition always starts 
from a decomposition of the target SIq- Nevertheless, in some special situa- 
tions it is still possible to achieve the patchy decomposition, using some a priori 
knowledge on the solution of the problem. This is the case of the infinite hori- 
zon problem associated to the linear- quadratic regulator, studied by Krener and 
Navasca in [3T] as a test for their patchy method (see also [30j and the Intro- 
duction for a brief description): 

min / (hyl + yl) + \aA dt subject to | f ^ f (13) 



aeR Jq \2 2 J [ J/2 = a 

with yi(t = 0) = xi and y2{t = 0) = X2. The exact value function for this 
problem turns out to be 



2\X2) \l V^, 

and it is easy to check that the origin (0, 0) is the only source of all characteristic 
curves. Then, using a small ball i?2(0,e) as a fictitious target, we are able to 
generate the patchy decomposition and then run the algorithm normally. In 
Fig. |6] we show the outcome of the simulation with^the following parameters: 
= [-1, l]2^ e = 0.05, i?=4, iVc=101, A = [-3, 3], iV=100, iV=200. Note that 
the choice R—A is due to limitations described in Section [4?5| Moreover, it is 
impossible to impose state constraints boundary conditions at the boundaries 



of the patches, since at some points the dynamics in (13) does not point inside 
the corresponding patch for any a G A. To fix this, wc imposed a Dirichlet 
boundary condition given by 

UpUv^U^p\dw. i = l,...,i?, (14) 

where Up^ is the first approximation of the solution computed in Step 2 of 
the patchy algorithm. Patchy errors are ||ii^p||i=0.002 and ||i?p||oo=0.009. Note 
that they are generally smaller than those in Tables l][3 (computed for other 



dynamics) because of the more favorable boundary condition ( 14 ) 



5 Patchy method's add-ons 

The patchy algorithm proposed in Section |3] has a multigrid nature, meaning 
that the computation of the solution on a rough grid is needed to start the 
optimal domain decomposition. Once this preliminary effort is done, it appears 



22 




(b) (c) 

Figure 6: (a) Decomposition in R—A patches, (b) error function and level sets 
of the patchy solution (the solution is truncated at value m = 3 in order to 
remove the boundary effects which are very important for this dynamics.), (c) 
patchy solution Up, (d) a detail of the error function shown in (b). 



to be natural to use all the information we have collected in order to speed up the 
proposed algorithm. First of all, we will always impose by default the boundary 
condition ( [l4| ) (note that it becomes available only after the computation on the 
rough grid). Further multigrid advantages we can take into account are listed 
in the following. 

AOl. We use Up^ computed in Step 2 as initial guess for Step 3.2. In this way 
we save some iterations to reach convergence. 

A02. Before Step 3.2 we order the nodes belonging to each patch in such a 
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way they fit as much as possible the causality principle |3S]. For exam- 
ple we can order the nodes with respect to their values. This ordering is 
optimal if the characteristic lines coincide with the gradient lines of the 
solution, as it happens in the case of the Eikonal equation. In general this 
is not true, anyway this ordering is often not too far from the optimal one. 



AOS. In Step 3.2 we reduce the number of discrete controls used in the numerical 
scheme, eliminating those controls which are "far" from the optimal one 
a*{xi) as computed by the first computation on the rough grid (Step 2). 
For example, if A — ^2(0, 1) we can introduce a reduction factor r > 1 
and replace A with the set 



a e A : a - a* > cos 



(7)} 



This is the only add-on which introduces a new error in the solution, which 
is anyway negligible in most cases. 

We point out that the patchy method can easily become an actual multigrid 
method. Indeed, we can in principle repeat the algorithm introducing a sequence 
of grids Gi, G2, . . . one finer than the other, until the desired precision is reached. 

In order to study the effect of the previously described add-ons, we introduce 
them one at a time and we compare the CPU time with the basic algorithm. 
Then, we apply all the features together. Results are reported in Table [9) 



Table 9: 
factor 4 



Effects of add-ons. 2 procs, R — 8, Nc — 32. Controls reduced by 



dynamics 


grid size 


basic 


AOl 


A02 


A03 


A01+A02+A03 


Eikonal 


100^ ^ 200^ 


20.0 


19.2 


9.6 


9.1 


5.7 


Eikonal 


100^ ^ 400^ 


130.7 


130.2 


40.5 


43.6 


17.8 


Eikonal 


100^ ^ 800^ 


928.1 


924.6 


238.8 


298.1 


100.6 


Fan 


100^ ^ 200^ 


31.9 


31.0 


11.4 


14.0 


7.6 


Fan 


100^ 400^ 


209.8 


205.7 


43.5 


72.3 


20.6 


Fan 


1002 ^ 8002 


1571.9 


1564.0 


247.3 


529.6 


110.6 


Zermclo 


100^ ^ 200^ 


23.2 


22.6 


11.5 


10.7 


6.7 


Zermclo 


100^ ^ 400^ 


143.5 


142.4 


46.2 


51.0 


20.3 


Zermelo 


100^ ^ 800^ 


1071.4 


1057.9 


290.1 


345.5 


111.3 



Note that CPU times for this test are lower than those in Sectionl4~3lbecause 



of the more favorable boundary condition ( 14 1 
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6 Numerical tests in dimension three 



We solve three 3D minimal time problems of the form ([ij. The numerical do- 
main is 51 = [—2, 2]-^ for all tests. 

Test 1 (Eikonal 3D) 

d = 3, f{xi,X2,X3,a)^a, A = 5,3(0,1) , flo = ^3(0, 0.5). 
Test 2 (Fan 3D) 

d^3, /(xi,a;2,a;3,a) = |a;i+a;2+a;3+0.1|a, ^ = ^3(0,1), llo = {a^i = 0}. 
Test 3 (Brockett problem |8ll29|) 

d = 3, f{xi,X2,X3,a) = {ai,a2,xia2-X2ai), A =[-5, 5]^, f^o = ^3(0, 0.25). 

For Tests 1 and 2 we used Nc— 189 discrete controls uniformly distributed on the 
unit sphere, while for Test 3 we used only Nc—9 discrete controls in {—5, 0, 5}^. 
The latter choice is motivated by the fact that using a larger number of discrete 
controls in [—5,5]^ do not lead to a different result. This is expected because 
here we solve the minimum time problem associated to the Brockett dynamics, 
then the optimal strategy always requires to saturate the control to the maximal 
admissible value (5, in this case). 

Results are reported in Table 10 Considering that for Tests 1-2 the number 



Table 10: 3D tests. 4 procs, R = 8. Add-ons enabled (Eikonal and Fan: controls 
reduced by factor 4. Brockett: not reduced) 



dynamics 


grid size 


CPU time 


II^^pIIi 


II^pIIoo 


Eikonal 3D 


50^ 100^ 


183 


0.033 


0.035 


Eikonal 3D 


50^ 200^ 


1217 


0.029 


0.042 


Fan 3D 


50^ ^ 100^ 


165 


0.064 


0.187 


Fan 3D 


50^ ^ 200^ 


1269 


0.056 


0.305 


Brockett 3D 


50^ -> 100^ 


132 


0.229 


0.024 


Brockett 3D 


50^ ^ 200^ 


1557 


0.165 


0.020 



of discrete controls is quite large, the CPU time is remarkable. Figure [7] shows 
the a level set of the value function for the Eikonal dynamics. It is perfectly 
visible the error located where the patches meet. Figure |8] shows the results for 
the Fan dynamics with 200^^ nodes. In Figs. [8]ja)-(b) we show the boundaries of 
the patches and some level sets of the solution, respectively. Level sets should 
be plans, but the state constraints imposed by the computational box Q bend 
them near 957. In Fig. l8](c) we show some optimal trajectories to the target. 
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Figure 7: One level set of the value function for the Eikonal dynamics 



Results for Brockett problem (Test 3) are different from the previous tests. 
First, CPU time turns out to be high with respect to the small number of dis- 
crete controls in use (just Nc — 9 controls). This could be related to the fact 
that characteristics are broken lines (see Figure [9](c)) that do not go directly 
to the target as in the Eikonal equation, nor bend slightly as for the Fan dy- 
namics (see Figure |8][c)), but change direction istantaneously (see also control 
switch regions in Figure [9](b)), so that this dynamics takes much more time 
to move information through the domain. As a consequence, it increases the 
number of iterations the scheme needs to reach convergence and also the time 
to compute the patches. On the other hand, the patchy error in norm is 
larger than the error in norm (see Tabic 10). This depends on the fact 
that the patchy decomposition obtained for this dynamics is rather complicate 
(see figure [9]( a)), in particular patches arrange themselves in (suggestive) sets 
with very large boundary areas, which increase the number of nodes where we 
commit the patchy error Ep. 
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